Constraining massive gravity with recent cosmological data 
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A covariant formulation of a theory with a massive graviton and no negative energy state has been 
recently proposed as an alternative to the usual General Relativity framework. For a spatially flat 
homogenous and isotropic universe, the theory introduces modified Friedmann equations where the 
standard matter term is supplemented by four effective fluids mimicking dust, cosmological constant, 
quintessence and stiff matter, respectively. We test the viability of this massive gravity formulation 
by contrasting its theoretical prediction to the Hubble diagram as traced by Type la Supernovae 
(SNela) and Gamma Ray Bursts (GRBs), the H{z) measurements from passively evolving galaxies, 
Baryon Acoustic Oscillations (BAOs) from galaxy surveys and the distance priors from the Cosmic 
Microwave Background Radiation (CMBR) anisotropy spectrum. It turns out that the model is 
indeed able to very well fit this large dataset thus offering a viable alternative to the usual dark 
energy framework. We finally set stringent constraints on its parameters also narrowing down the 
allowed range for the graviton mass. 
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I. INTRODUCTION 



There is a limited but crucial number of cosmological 
evidences which does not fit in the scheme constituted by 
General Relativity (GR) and the Standard Model of par- 
ticle physics. On the one hand, it seems that the amount 
of baryonic matter in galaxies and galaxy clusters is not 
sufficient to cause the observed behaviour of the gravita- 
tional field on those scales. On the other hand, recent 
data on the SNela Hubble diagram provide evidences 
for an accelerated expansion, also confirmed by several 
other independent cosmological observations [HHI . Such 
a large dataset can be very well reproduced by the con- 
cordance ACDM model made out of cold dark matter 
(CDM) and cosmological constant A. Notwithstanding 
this remarkable success, the ACDM scenario is theoreti- 
cally unappealing since it comprises two exotic matter 
sources, namely dark energy and dark matter, and is 
moreover plagued by several well known shortcomings. 

A different constructive point of view is to address 
these problems introducing modifications of GR over 
large distances in order to possibly explain our ignorance 
about the 95% of the energy and matter content of the 
Universe. In this setting, a massive deformation of GR is 
a plausible modified theory of gravity that is both phc- 
nomenologically and theoretically intriguing. 
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From the theoretical point of view, a small nonvanish- 
ing graviton mass is an open issue. The idea was origi- 
nally introduced in the work of Fierz and Pauli 0] , who 
constructed a massive theory of gravity in a flat back- 
ground that is ghost - free at the linearized level. Since 
then, a great effort has been put in extending the result at 
the nonlinear level and constructing a consistent theory. 
As a somewhat unanticipated result, we nowadays recog- 
nize that such models can be relevant for cosmology too 
as possible alternative candidates to drive the accelerated 
expansion without the need of any exotic component. 

Here we study cosmological solutions in the frame- 
work of a covariant massive gravity model, recently pro- 
posed in [1, [f|. At the linearized order, the mass term 
breaks the gauge invariance of GR. Moreover, in order 
to construct a consistent theory, nonlinear terms should 
be tuned to remove order by order the negative energy 
state in the spectrum Q. The theoretical model under 
investigation follows from a procedure originally outlined 
in 0, H| and has been found not to show ghosts at least 
up to quartic order in the nonlincarities [9|, [1^ . 

We focus on cosmological equations of such theory, 
considering a spatially flat Robertson- Walker (RW) 
model in presence of both matter and radiation. Due 
to the potential term, four new terms are present in the 
Friedmann equations mimicking dust, cosmological con- 
stant, quintessence and stiff matter. The presence of the 
cosmological constant and quintessence - like terms easily 
suggest the possibility to achieve an accelerated expan- 
sion, while the stiff matter one can guarantee the halting 
of the cosmic speed up in the intermediate redshift regime 
thus recovering the standard decelerating matter domi- 
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nated epoch. In order to check whether this is indeed 
the case, we contrast the model predictions with a wide 
dataset thus also being able to constrain its parameters. 
As a side effect, this will also allow us to narrow down 
the range for the graviton mass. 

The plan of the paper is as follows. Sect. II presents 
the basics of the theory and put down the modified Fried- 
mann equations describing the background cosmic evo- 
lution. The data we use and some subtleties specific to 
the massive gravity framework are discussed in Sect. Ill, 
while the results of the likelihood analysis and the con- 
straints on the massive gravity coupling quantities are 
given in Sect. IV. We finally summarize and highlight 
some important issues in the concluding Sect. V. 



@. For generic values of the (03,04) couplings and a 
spatially flat universe consistent with CMBR data fl2l ]. 
the expansion rate turns out to be (lH ] 



3H 2 
8^G 



= PAi(a) + p r (a) 

+ m 2 {4c 3 + c 4 - 6) 
3m 2 /3(3 - 3c 3 - c 4 ) 



+ 



3m 2 /3 2 (c 4 + 2c3 - I) 
m 2 /3 3 (c 3 + c 4 ) 



(3) 



II. MASSIVE GRAVITY 

Let us consider a four dimensional manifold equipped 
with both a dynamical metric g M „ and a non - dynamical 
flat metric f^. In order to determine the dynamics on 
this structure, the starting point is to consider the fol- 
lowing action 1 



S = -8^J{l R + m2U ) d4x + S " 



(1) 



where G is the Newton coupling constant, R the Ricci 
scalar for g^ v , IA the potential, and Sm describes ordinary 
matter which is supposed to directly interact only with 
g^ v . The potential term, coupled through the graviton 
mass to, is given by 



1 

2 



(2) 



where a is the scale factor, H — a/a the usual Hubble pa- 
rameter (with the dot denoting derivative with respect to 
cosmic time), to 2 = (m 2 c 2 )/(8irG) and we have assumed 
a source term made out by dust (M) and radiation (r). 
Since matter is assumed to be minimally coupled to grav- 
ity only, its conservation equation will be the usual 



3H( Pi + Pi )=0 



(4) 



with pi = (pi = p r /3) for i = M (i = r). Introducing 
the redshift z = 1/a — 1 , we can therefore conveniently 
rewrite Eq.(j3]) as 



eff 



(i 



(5) 



with E(z) = H(z)/Hq (the label denoting present day 
values), fli the standard present day density parameter 
for the component i, while we have defined the effective 
density parameters for the massive gravity terms as 



with 



the Levi-Civita symbol, while the field K, is 



related to the flat metric through 



and we have defined 



7^ 



The theory is fully assigned by the graviton mass m g = 
hm/c and the two coupling parameters (03, c 4 ). 

In order to derive the cosmological equations, one has 
to insert the RW metric into the field equations obtained 
by the usual variational approach to the action (TTJ and 
impose the constraints dictated by the Bianchi identities 



f n e A ff =/i 2 (4c 3 + c 4 -6) 



n e ^ f = 3/i 2 /3(3 - 3c 3 - c 4 ) 
Q e s ff = 3fl 2 (3 2 (c 4 + 2c 3 - 1) 

I n% f = n M -m 2 /3 3 (c 3 + Ci ) 



(G) 



with fx = m / Pcrit and p c 



3H 2 f 



inG the present 



day critical density 2 . Note that the (c 3 , c 4 , fi) parameters 
entering Eq.Q are all dimcnsionless. 

A look at Eq.([5]) shows that a massive graviton has a 
double effect on the Hubble rate. First, as is somewhat 
expected, it changes the matter content shifting the dust 



1 Unless otherwise specified, we will set the speed of light c = 1. 



2 Note that, since in = cm g /h, it has the dimension of (length) 1 
so that p, 2 = (m 2 c 2 )/ (3H 2 ) is dimcnsionless. 



3 



density parameter from the usual s!m value to the effec- 
tive one £1%^ ■ It is, however, worth stressing that, since 
the sign of C3 + C4 is not set from the theory, it is also 
possible that the effective matter content is smaller than 
the actual one. Actually, the most interesting feature is 
the presence of the three further terms (Q e /^ , flg^ ,fl e /^) 
mimicking a cosmological constant, a quintessence - like 
field and stiff matter, respectively. Depending on the val- 
ues of (03, C4) and the graviton mass, it is therefore pos- 
sible to reproduce the expected series of radiation domi- 
nation, matter era and accelerating expansion (driven by 
both and fig^) with only a small contribute from 

the unusual stiff matter term. It is this particular feature 
that makes massive gravity so appealing from a cosmo- 
logical point of view motivating the present analysis. 

In order to constrain the model parameters, it is actu- 
ally more convenient to replace the effective density pa- 
rameters with wo = w e ff(z = 0) and w p = dw e ff/dz\ z= o 
with w e ff the equation of state (EoS) of the effective 
dark energy model having the same Hubble parameter 
than the massive gravity one. This quantity is defined as 



W eff 0) = -1 + 



2 d\nE(z) 
3dln(l + z) 



n e //{i + z) 3 n r (i + z) 4 

E 2 (z) E 2 {z) 

n e //{i + z ) 3 n r (i + z) 4 



E 2 (z) 



so that, by using Eq.©, we get 



E 2 {z) 



(7) 



(/>,, 



ff (z) = -{s[n e / f + n e Q "{i + z) + nf f (i + z f]} 



zn e / f + 20^(1 + z) 



(8) 



Imposing now E 2 (z = 0) = 1, Wo = w e ff(z — 0) and 
w p = dw e f f I 'dz\ z= Q, it is only a matter of algebra to get 



O 



eff 



9, 










— QSI^Wq 




- fl eff - 


- Q r )w p 


+ 3(1- 


- n eff ~ 


2Sl r ) , 



eff 



1 - Q 



eff 



eff 



eff 
M 



(10) 



(11) 



Using these relations, we can conveniently parameterize 
the model in terms of the effective matter density pa- 
rameter fi^ and the present day values of the effective 
EoS and its derivative (wq,w p ). It is worth stressing 
that, depending on (fl e J/ , f2 7 , wo, w p ) combinations, it is 
possible that some of the (fi^^, fig , fig^) parameters 
take a negative value. This is, however, not a problem 
since (n^ , f2g^, ) refer to effective fluids, not ac- 
tual ones. This can also be seen from Eqs.(|6]) which show 
that negative fi^^ can indeed be obtained depending on 
the values of the massive gravity couplings (03,04). 



III. MASSIVE GRAVITY VS DATA 

Eq.© shows that the massive gravity introduces an 
effective dark energy fluid whose EoS has the right be- 
haviour to drive cosmic acceleration. Indeed, should the 
f^A t erm be the dominant at low z, we get w e ff(z << 
1) ~ —1 so that a ACDM-like expansion is obtained in 
the late universe. As z increases, we can still have an 
accelerated expansion as far as £Iq^ dominates leading 
to w e ff(z) ~ —2/3 in the intermediate redshift regime, 
while the transition to a decelerating epoch is guaran- 
teed by the stiff matter term f2g^. These encouraging 
features have to substantiated by the comparison with 
the available data which can both validate the model and 
constrain its parameters. 

Although the dataset we are going to use is a typical 
one, there are some subtleties specific of the massive grav- 
ity framework which require some caution. We therefore 
prefer to spend some words to describe how we compare 
the model to the each kind of dataset. 



A. Hubble diagram 

= 9(1 — ^a/ ~ ^r)wo A first step in testing any proposed model is the com- 

parison with the Hubble diagram, i.e. the distance mod- 
ulus fj, as a function of the redshift z. This is related to 
the Hubble parameter E(z) as 

^ = 25 + 51ogd L (z,p) (12) 



-!<> 


- n e // - n r )w 2 Q 






- 'iQl 1 / - AQ r )w 




+ 2* 1 






+ (1- 




(9) 



4 



with 



d L {z,p) = (1 + z)r(z,p) 
the luminosity distance and 



(13) 



C 



GRB 



exp 



where 



X 2 GRB (P) 



(19) 



r(z,p) 



dz' 



Ho J E(z',p) 



(14) 



the comoving distance. Note that we have collectively 
denoted with p the set of model parameters. 

A classical tracer of the Hubble diagram is represented 
by SNela so that we use the Union2 dataset 14| compris- 
ing AfsNeia = 557 objects tracing the Hubble diagram 
over the rcdshift range (0.015, 1.4). The corresponding 
likelihood term will read : 



1 



-SNela 



(2tt)Mw.i.|C 



SNela 



11/2 



(15) 



x exp 



^SNelaiP^SNela^SNelaip) 



with C sNeia the SNela covariance matrix and D(p) a 
NsNeia - dimensional vector with the i - th element given 
by /J-obs(zi) — ^th{zi, p), i.e. the difference between the ob- 
served and predicted distance modulus. Neglecting sys- 
tcmatics, the covariance matrix becomes diagonal and we 
can simplify Eq. (fT5|) as 



C s 



1 



Nela 



1/2 

SNela 



exp 



XsNelafa) 



(16) 



where 



XsNela 



E 



Hobs(Zi) - Hth{Zi, p) 



FsNela = 

i=l 



(17) 



(18) 



While SNela efficiently probe the late universe thus being 
highly sensible to the present day cosmic speed up, the 
transition to the decelerated expansion and the matter 
dominated era are better investigated resorting to higher 
z tracer. This is provided by GRBs so that we rely on 
their Hubble diagram as derived in [l5| based on a model 
independent calibration of five different scaling relations 
and the catalog given in [l6[ . Note that we cut the sample 
only using Mgrb = 64 objects probing the range 1.48 < 
z < 5.60 in order to avoid any residual correlations with 
the SNela sample (which is used to calibrate the GRBs 
scaling relations) . The GRB likelihood term is similar to 
the SNela being defined as 



Xgrb 



Ngrb 

E 



GRB 



^obs(z l ) - pL t h{Zi,p) 



Mgrb 

n 

i=l 



a int) 



(20) 



(21) 



Note that we have here added to the measurement un- 
certainties <Tj the intrinsic scatter Oi n t which takes into 
account the dispersion of the single GRBs around the in- 
put scaling relations. We will marginalize over Oi nt so 
that we will not give any constraint on this quantity. 



B. H(z) data 

Being a probe of the integrated Hubble parameter, the 
SNela + GRB Hubble diagram smooths out the devia- 
tions from the (unknown) best fit model provided these 
are not too extreme. It is therefore desirable to add a 
second dataset which directly probes the E{z) which is 
indeed possible resorting to the differential age method 
fl7| . Such a technique is motivated by noting that, from 
the relation dt/dz = —(1 + z)H{z), a measurement of 
dt/dz at different z gives, in principle, the Hubble param- 
eter. Good results can be obtained by using fair samples 
of passively evolving galaxies with similar metallicity and 
low star formation rate so that they can be taken to be 
the oldest objects at a given z. Stern et al. [l8| have 
then used red envelope galaxies as cosmic chronometers 
determining their ages from high quality Keck spectra 
and applied the differential age method to estimate H (z) 
over the redshift range 0.10 < z < 1.75 [nj. We use their 
data as input to the following likelihood function : 



Ch(p) 



(27r)^/ 2 r^ /2 



exp 



X 2 H (P) 



(22) 



with 



* 2 * = E 



H bs{z.i) - H(zi,p) 



i=l 

where Mh = 11 is the number of sample points. 



(23) 



(24) 
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C. Baryon Acoustic Oscillations 



An alternative and very promising method to trace 
the distance - redshift relation relies on the measurement 
of baryon acoustic oscillations (BAOs) in the large scale 
clustering pattern of galaxies [20(. BAOs correspond to 
a preferred length scale imprinted in the distribution of 
photons and baryons by the propagation of sound waves 
in the relativistic plasma of the early universe. This 
length scale, corresponding to the sound horizon r s (zd) 
at the baryon drag epoch, manifests itself in the cluster- 
ing pattern of galaxies as a small preference for pairs of 
galaxies to be separated by r s (zd), causing a distinctive 
peak in the 2 - point galaxy correlation function. Since 
the Fourier transform of a 5 - function is a sine, this sig- 
nature will look like a series of decaying oscillations in the 
galaxy power spectrum thus motivating the BAO name. 
The small amplitude of the peak and the large size of 
the relevant scales imply that large volumes (~ 1 Gpc 3 ) 
and a high number of galaxies (~ 10 5 ) must be observed 
in order to ensure a robust detection. It is therefore not 
surprising that the first reliable determination plj has to 
await for the third data release of the SDSS survey with 
other determinations f22j all relying on similarly large 
(both spectroscopic and photometric) galaxy surveys. 

In order to include BAOs as constraint, we use 
the measurement at six different redshifts, namely 
z = 0.106 from 6dFGRS [H, z = (0.20,0.35) from 
SDSS + 2dFGRS 0, and z = (0.44,0.60,0.73) from 
WiggleZ J25(. Since the three surveys are independent 
on each other, we can define the likelihood function as 



£-BAO = £edFGRS x &SDSS x ^-WiggleZ 

with the three terms given by [25| 



6dFGRS 



v/2^2 



exp 



!()(> 



^0.106 ^0.106 (P) 



CO. 106 



c 



1 



: (2tt)"soss\ Csdss \^ 

T) SDSS(P) C SDSS' D SDSS(P) 



(25) 



(26) 



x exp 



, r B (zd) , , 
d z = , , . = r s {z d ) x 

dy \Z) 



czr 2 [z) 



H E(z) 



(28) 



with the sound horizon to distance z given by 



r s{z) 



E- 1 (z)dz' 



o Jz v/l + (3w b )/(4w r )(l + z')- 



(29) 



Note that we have here introduced the physical baryon 
and photon density parameters defined as Ui = f2j/i 2 
(with h the Hubble constant in units of 100 km/s/M pc) , 
while Zd is the drag epoch redshift. We will follow [23[ 
to set (dg h f 06 , ao.io6)j while [24| gives the relevant d° bs 
values and covariance matrix for the measurement at 
z = (0.20, 0.35) from SDSS + 2dFGRS data. 

The WiggleZ survey has recommended to use a differ- 
ent BAO related quantity so that DwiggieZ is a three 
dimensional vector whose i - th element is given by the 
difference between the observed and predicted value of 
the acoustic parameter A(z) defined as [2l|, [25[ 



A(z) 



y/fl M H§d v (z) 



(30) 



with the volume distance dy{z) given in Eg. (|2"5)) . We use 
the observed values and their covariance matrix for A(z) 
determinations at z = (0.44,0.60,0.70) reported in |25j . 

An important remark is in order here. Both the drag 
epoch redshift z d and the acoustic parameter A(z) ex- 
plicitly depend on the matter physical density parameter 
u>M- It is worth wondering which value has to be used 
for Qm in the definition of the i.e. whether it is 
u>m = ilj\f/i 2 or ojm = il e // h 2 . This is a subtle question 
that has not a definitive answer at the moment. Indeed, 
both Zd and A(z) arc related to the clustering properties 
of the matter so that one should investigate how clus- 
tering takes place in massive gravity to understand if it 
depends on the amount of either effective or actual mat- 
ter. In order to be conservative, we have decided to set 
lu]\i = Q,Mh 2 when computing both Zd and A(z) since 
the standard dust matter is the only one which surely 
clusters. Actually, since the graviton mass is expected to 
be very small, we can anticipate that £Im and Q e // will 
not be very different so that, should our choice be the 
incorrect one, the estimated z d and A{z) values will be 
biased low by a very small amount. 



C 



1 



WiggleZ 



(27r) A/ ""3s i = z |C 



WiggleZ 



11/2 



^WiggleZ (V^WiggleZ^* 



(27) 

WiggleZ {p) 



exp 



In Eq. (|27p . D$dss is a two dimensional vector with the 
values of d° bs — d t J l (p) having defined 



D. CMBR data 

Combining the data described above, we are able to 
trace the background expansion of the late and interme- 
diate z universe. CMBR data, on the contrary, give us 
a picture of the universe in its infancy, at the redshift of 
the last scattering surface (z* ~ 1100). As shown in [l2j], 
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rather than fitting the full CMBR anisotropy spectrum, 
one can get reliable and almost equivalent constraints by 
considering the so called distance priors, namely the red- 
shift z* of the last scattering surface, the acoustic scale 
I A = nr(z+)/r s (z+) and the shift parameter 1Z [2"o| 



c/H 



(31) 



The CMBR likelihood function is then defined as 



C 



1 



CMBR 



(27r)^MB R |c CMBfl |V2 

■ 

x exp 



(32) 



Dcmbr(p)CcmbbJ-* 



CMBR 



(P) 



where the Ncmbr{= 3)-dimensional vector Dqmbr con- 
tains the difference among observed and theoretically 
predicted values of (£a,7Z, z±). We take the observed 
distance priors and the corresponding covariance matrix 
from [l2j and use the approximated formula in (2?J to 
compute z+ as function of (ui , uj^)- We stress that we set 
here ujm = ^Mh 2 , but use £l e ^ in the TZ estimate. This 
latter choice has been motivated by the consideration 
that the shift parameter actually turns out by neglecting 
terms other than the (1 + z) 3 in the evaluation of the 
distance to the last scattering surface. Since the impor- 
tance of this term is determined by the f2j$ parameter, 
it is this quantity which must enter the TZ definition. 



E. The full likelihood function 

Since the single datasets are independent on each 
other, combining all of them in a single fit is straight- 
forward. One has only to define a full likelihood as 



(fl e // ,wo,Wp, H 7 ujm), but leading to unphysical massive 
gravity parameters. Indeed, while Eqs.([§| - (fTTj) always 
give real values whatever the (fi^^, Wq, w p ) parameters 
are, solving Eqs.dB]) with respect to (p, 2 , /?, C3, C4) can lead 
to unphysical values such as a negative to 2 . In order to 
avoid this possibility, we therefore define our final likeli- 
hood function as 



£(P) = Cdatafa) X V(p) 



(33) 



with V(p) = 1 if i.) p? > 0, m.)/3 > (in order to avoid a 
change in the signature of the metric) and Hi.) C3C4 < 
[29l ] . while it is 'P(p) = otherwise. 

In order to efficiently explore the parameter space, we 
use a Markov Chain Monte Carlo (MCMC) method run- 
ning three parallel chains and keep adding points until 
the Gelman - Rubin [30j convergence criterium is satis- 
fied. The best fit model will be the one maximizing the 
likelihood function £(p), but, in a Bayesian framework, 
most reliable constraints on a single parameter pj are 
given by the mean, median and 68% and 95% confidence 
levels (CL) estimated from the likelihood function after 
marginalizing over all the parameters but the j - th one. 
Actually, the marginalization is implemented by simply 
looking at the histograms of the merged chain after cut- 
ting the burn - in phase and thinning to avoid spurious 
correlations. 



IV. RESULTS 

Three parallel chains with ~ 13000 points each turn 
out to be sufficient to achieve a good convergence of the 
MCMC code, the R parameter being smaller than 1.1 for 
all the fitted parameters. The best fit model read 



n eff 



0.282 , ui M = 0.1351 



0.697 



C-data — C-SNela x C-GRB x C-H x C-BAO x &CMBR x A] 

where the last term is a Gaussian prior on h given by 
the SHOES collaboration [28| which have determined 
h = 0.738 ± 0.024 from the local distance ladder. 

According to the above review of the data, Cdata will 
depend on the massive gravity parameters (ft e J/ ,wq, w p ), 
the present day Hubble constant h and the physical den- 
sity parameters (cob, ujm, oj t ) of baryons, matter and ra- 
diation. In order to reduce the parameter space, we will 
set (w6 i w 7 ) = (2.258 x 10" 2 , 2.469 x 10~ 5 ) in agreement 
with [12J and set the radiation density parameter as 

n r = cj 7 /i- 2 (1 + 0.22717V e// ) 

with N e ff = 3.04 the effective neutrino number. 

Should we rely on Cdata only, nothing would prevent 
the likelihood to favour models with reasonable values of 



w = -1-04 , w p = 0.00 , 

while its theoretical predictions for the Hubble diagram 
and the Hubble parameter are superimposed to the data 
in Fig.[TJ Here, we also plot the effective dark energy 
EoS for the best fit model and the constraints from the 
MCMC analysis. To this end, for each z, we evaluate 
w e ff(z) for all the points in the chain and infer the me- 
dian value and 68% confidence ranges from the result- 
ing distribution. Note that the best fit w e //(z) and the 
median thus obtained are not equal because of the de- 
generacies in the model parameters space so that it is 
actually the median w e ff(z) the most reliable one in a 
Bayesian framework since it takes explicitly into account 
such degeneracies and marginalize over them. 

It is interesting to look in detail to how the best fit 
model compares to the observed data to guess some hint 
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FIG. 1: Left. Best fit model superimposed to the Hubble diagram. Centre. Best fit model overlapped with the H(z) data. 
Right. Effective EoS with the solid line referring to the best fit model and the dashed ones to the median and 68% CL. 



Id xbf (x) x 68% CL 95% CL 

Q.\{ s 0.282 0.284 0.284 (0.272, 0.293) (0.260, 0.306) 

lum 0.1351 0.1375 0.1380 (0.1331, 0.1412) (0.1274, 0.1479) 

h 0.697 0.705 0.703 (0.695, 0.714) (0.689, 0.725) 

wo -1-04 -1.09 -1.09 (-1.13, -1.04) (-1.17, -0.98) 



is parameterized by (wo,w p ) = (—1,0). The smoothing 
process introduced by the integration entering the lumi- 
nosity distance definition cancels out the residual differ- 
ences which, on the contrary, are better appreciated when 
fitting the H(z) data leading to d.o.f. = 1.23. This 
is still a fully acceptable result considering both the low 
statistics (giving d.o.f. = 5) and the scatter of the H{z) 
data around the best fit ACDM prediction. 

The model turns out to be well efficient in reproducing 
also the BAO data. For the distance ratios, we get 



w p 0.00 -0.05 -0.05 (-0.09, -0.01) (-0.11, 0.02) 



TABLE I: Constraints on the fit parameters. Columns are as 
follows : 1.) id parameter, 2.) best fit, 3.) mean, 4.) median, 
5.), 6.) 68% and 95% confidence ranges. 



d b f 1Q6 = 0.3435 vs d%% 6 = 0.336 ±0.015 
4/.200 = °- 1872 vs d o%a = °- 1905 ± °- 0061 > 



^350 = 0.1123 vs d%% = 0.1097 ±0.0036 



on pro and cons of the massive gravity framework. In- 
deed, while the low reduced x 2 values make us forecast 
that the agreement will be very good, it is nevertheless 
worth checking what happens for each single probe to 
see whether the likelihood is dominated by a single term 
only 3 . First, we note that the model very well fits the 
Hubble diagram data with 

XsNeiJd-o.f. = 0.98 , X grb / d.o.f . = IM , 

d.o.f. = Hi — n p being the number of degrees of free- 
dom. This is a nice consequence of the effective EoS 
turning out to be very close to the ACDM one, which 



3 To understand why this could happen, one can imagine to fit 
only the SNela and GRBs Hubble diagrams. Since there are 
an order of magnitude more SNela than GRBs, should the fit 
to the two separated datasets points towards different regions 
of the parameter space, the fit to the combined sample will be 
closer to the SNela only one. In such a case, one could still find 
a low reduced Xgrb because of the larger uncertainties on this 
dataset. Actually, this does not take place, but we have preferred 
to explicitly check for this possibility. 



while the acoustic parameter values read 



' „4 b /(0.44) = 0.460 vs .A o6s (0.44) = 0.474 ±0.034 
< .4^ (0.60) = 0.435 vs _4 obs (0.60) = 0.442 ±0.020 . 
t A bf (0.73) = 0.414 vs A obs (0.73) = 0.424 ±0.021 

As it is apparent, the theoretically predicted values agree 
with the observed one well within lcr so that we can be 
well confident in the ability of the model to give the cor- 
rect BAO features. It is worth noting, however, that 
the acoustic parameter values are systematically smaller 
than the observed one with A bf (z)/A obs (z) ~ 0.98 in- 
dependent on z. Noting that, for the best fit values, 
= 0.99, one can argue that the observed dis- 
crepancy could be remedied by assuming that the full 
effective matter term enters the clustering process or, at 
least, should be taken into account when computing the 
BAO acoustic parameter at different redshifts. 

Finally, we compare the theoretically predicted dis- 
tance priors with the observed ones being 
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( i h [ = 302.023 vs i°^ s = 302.09 ± 0.76 



K bf = 1.729 vs TZ obs = 1.725 ± 0.018 



z h J = 1090.92 



vs zf s = 1091.30 ±0.91 



so that the agreement is still excellent. 

The above results have been obtained contrasting the 
model to the full observational probes we have described 
in Sect. Ill which also include data from GRBs and the 
H (z) as inferred from galaxy ages. Although both these 
datasets are well fitted by the model, their use as cosmo- 
logical probes is quite recent so that one can not exclude 
a priori that some undetected systematic is present and 
potentially bias the constraints. As a consistency check, 
we have therefore repeated the above analysis removing 
both GRBs and H(z) data from the fit. The median 
value and 68% CL turn out to be : 



o e ff — n oqo+0.012 

si M — u - zoz -0.012 



1316+ 0039 
u.ioio_ a0061 



0.685 



-0.012 



w a = -1.061°;°;? 



= 0.06 



+0.03 



These results are in good agreement with those in TableU] 
with the 68% CL well overlapping for all the parameters 
but w p . This latter quantity is the one with the largest 
shit since now positive rather than negative values are 
preferred. Actually, for both fits, the results argue in 
favour of a very mild evolution of the effective EoS which 
can be easily smoothed out from the integration needed 
to get the luminosity distance. We therefore expect that 
adding H(z) data erases this discrepancy. Indeed, fitting 
this dataset too (but still excluding GRBs), we get : 



ft e i f = 0.284t - 008 



71 n+ 005 
u. i iu_ 010 



-0.006 



lum = 0.1399 



+0.0023 
-0.0061 ) 



wo = -1.07^ 



+0.03 
05 



-o r)7 +a03 

u - u ' -0.03 



which are now in still better agreement with those in Ta- 
bleU] with the w p confidence range moving towards neg- 
ative values. We also note that the median h value is 
larger and almost equal to the one from the fit to the full 
dataset. This comes out as a consequence of the need 
for the fit to align the present day Hq value to the trend 
followed by the z > measurements. That this is the 
case is also confirmed by the last test we do removing 
the H(z) data, but adding the GRB ones. We now get : 



M 



0.28T 



-0.011 
-0.012 



uj m = 0.1347; 



0.695 



-0.010 
-0.005 i 



w 



1.07 



-0.04 
-0.06 



-0.0048 
-0.0055 ' 



-o.or 



0.03 



Id xbf (x) x 68% CL 95% CL 

m g 2.85 7.26 0.06 (0.02, 0.96) (0.01, 30.6) 

13 0.39 0.38 0.41 (0.16, 0.54) (0.05, 0.62) 

c 3 2.42 5.91 3.86 (1.57, 11.2) (1.31, 21.2) 

c 4 -2.88 -8.15 -5.52 (-15.8, -0.82) (-29.8, -0.70) 



TABLE II: Constraints on the massive gravity parameters 
with m g the graviton mass in 10 -66 g units. 



We indeed find a smaller h value, but still a good agree- 
ment with the results in TableU] We therefore conclude 
that, should some undetected systematics be present in 
the GRBs and/or H(z) data, the constraints would not 
be affected so that we can safely rely on them. 

Having convincingly shown that the model is able to 
reproduce a wide set of data, it is worth looking at the 
constraints on the quantities explicitly entering the mas- 
sive gravity Lagrangian. For each set of fitted parameters 
along the merged MCMC chain, we therefore solve for the 
(/i 2 , /3, C3, C4) and give the inferred constraints in TablelTTl 
where we report the graviton mass m g = 2>Hoflh/c in- 
stead of jl since the former is more interesting from a 
theoretical point of view. 

A caveat is in order here. Solving for (/i 2 , /3, C3, C4) as 
function of the fitted parameters can give rise to large val- 
ues of the (03, C4) couplings for particular combinations of 
the input quantities. As a consequence, the (C3, C4) distri- 
butions turn out to be strongly asymmetric with ~ 90% 
of the points being close to the median value, but the 
remaining ~ 10% distributed in a very long tail towards 
extremal values (positive for C3 and negative for C4). A 
similar argument also applies for m g which has a long 
tail towards large masses. In order to avoid this effect, 
we have arbitrarily cut these tails rejecting the most ex- 
treme 10% points and recomputing the constraints from 
the survival points. These are the values reported in Ta- 
ble|TT] which can be considered as a compromise between 
what data tell us and what is actually reasonable. It is 
worth stressing that the above cut has only a minor effect 
on the constraints on the graviton mass. 



V. CONCLUSIONS 

Until not long ago, it was thought to be impossible 
to construct a consistent theory of massive gravity com- 
patible with current observations [1, [3l|. Only in 2010, 
de Rham and Gabadadze [1] presented a good candidate 
for a ghost - free completion of the Fierz - Pauli theory by 
choosing interactions with no more than two time deriva- 
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tives of the scalar degrees of freedom [4J . The cosmologi- 
cal sector of this theory have been recently worked out in 
0, [HI thus setting the necessary framework to compare 
the model predictions with the observed universe. We 
have then been able to constrain the model parameters 
using a large dataset comprising the SNcIa + GRB Hub- 
ble diagram, H(z) measurements from cosmic chronome- 
ters, BAOs data and the CMBR distance priors. It turns 
out that the model is in very good agreement with the 
data thus giving an observationally motivated support to 
massive gravity as a theoretically appealing alternative 
the the troublesome ACDM scenario. 

Lacking up to now a clear physical interpretation of 
the (c3,Ci) couplings entering the modified gravity La- 
grangian, one can only indirectly check which are the 
conditions they have to fulfil in order the theory be con- 
sistent. To this end, we stress that our constraints are 
consistent with the conditions C3 + C4 = and C4 7^ — f , 
but the case C3 + C4 > is still allowed. In both cases, 
the constrained parameter space refers to models which 
have the right properties to recover GR in the low energy 
limit as demonstrated by [34[ through the study of spher- 
ically symmetric solutions. We therefore end up with a 
massive gravity theory able to reproduce the cosmolog- 
ical data without altering the success of the standard 
General Relativity on the Solar System scales. 

As a further relevant result, the comparison of the the- 
ory with the cosmological data has allowed us to set a 
limit on the graviton mass. Although an explicit deriva- 
tion is needed, we can assume that, in the low energy 
limit, the gravitational potential will be boosted by a 



Yukawa -like correction modulated by the length scale 
X g oc \/m g . One can then compare our results with 
the upper limit m g < 7.68 x 10~ 55 g from the dynamics 
in the Solar System (32j and the more stringent limit, 
m g < 10 -59 g, derived by requiring the derived dynam- 
ical properties of a galactic disk to be consistent with 
observations [33| . Our 95% upper limit is six orders of 
magnitude smaller than the latter one. We can therefore 
conclude that the agreement with the cosmological data 
leads to a massive gravity theory which has no impact 
on the dynamics of sytems in the low energy regime thus 
avoiding any difficulty in matching different scales. 

Finally, since covariant massive gravity seems to be a 
consistent modification of GR, it would be interesting and 
worthwhile to further explore it. In particular, investiga- 
tion should be improved by analyzing perturbed solutions 
and full prediction of cosmological perturbations. 
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